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Abstract: We study the generation of a kination-dominated phase in the context of a 
quintessential model with an inverse-power-law potential and a Hubble-induced mass 
Q-V term for the quintessence field. The presence of kination is associated with an oscillating 

evolution of the quintessence field and the barotropic index. We find that, in sizeable 
regions of the parameter space, a tracker scaling solution can be reached sufficiently 
early to alleviate the coincidence problem. Other observational constraints originating 
£f) ■ from nucleosynthesis, the inflationary scale, the present acceleration of the universe 

and the dark-energy-density parameter can be also met. The impact of this modified 



> 



t^J- \ kination-dominated phase on the thermal abundance of cold dark matter candidates is 

investigated too. We find that: (i) the enhancement of the relic abundance of the WIMPs 
with respect to the standard paradigm, crucially depends on the hierarchy between the 
freeze-out temperature and the temperature at which the extrema in the evolution of the 
quintessence field are encountered, and (ii) the relic abundance of e-WIMPs takes its 
present value close to the temperature at which the earliest extremum of the evolution 



of the quintessence field occurs and, as a consequence, both gravitinos and axinos arise 
as natural cold dark matter candidates. In the case of unstable gravitinos, the gravitino 
constraint can be satisfied for values of the initial temperature well above those required 
in the standard cosmology. 
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1. Introduction 

A plethora of recent data [1,2] indicates [3] that the two major components of the present universe 
are Cold Dark Matter (CDM) and Dark Energy (DE) with density parameters [1]: 

(a) Ocdm = 0.214 ± 0.027 and (b) Q D e = 0.742 ± 0.03, (1.1) 

at 95% confidence level (c.L). Identifying the nature of these two unknown substances is one of the 
major challenges in modern cosmo-particle theories. 

Among the natural candidates [4] to account for CDM are [5] the weakly interacting massive 
particles (WIMPs) [6, 7], with prominent representative in supersymmetric (SUSY) theories the 
lightest neutralino [8] and the extremely WIMPs (e-WIMPs) [9] with most popular representatives, 
the gravitino, G, and the axino, a. Assuming /^-parity conservation, the lightest SUSY particle 
(LSP) is stable and can be either a WIMP or an e-WIMP, in a sizeable region of the SUSY pa- 
rameter space. The interactions of WIMPs ensure that they come to chemical equilibrium with 
the plasma and decouple from it at a temperature 7> ~ (10 - 20) GeV. On the other hand, the 
interaction of e-WIMPs are extremely weak since they are suppressed by the reduced Planck scale, 
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m P = M P / V8tt (where M P = 1.22 • 10 19 GeV is the Planck mass) in the case of G and by the 
axion decay constant, f a ~ (10 10 - 10 12 ) GeV (for a review, see Ref. [10]) in the case of a. Con- 
sequently, e-WIMPs depart from chemical equilibrium very early and their relic density (created 
due to this early decoupling) is diluted by primordial inflation. Subsequently, e-WIMPS can be re- 
produced in the thermal bath through scatterings [11, 12, 13, 14] and decays [13, 15, 16] involving 
superpartners. For both WIMPs and e-WIMPs (hereafter denoted collectively as X) we may have 
extra non-thermal contributions to their relic density, Clxh 2 , [17, 18], from the out-of-equilibrium 
decay of the next-to-LSP (NLSP). However, this mechanism is highly model dependent since it is 
sensitive to the nature of the NLSP, its decay products and the extra restrictions which have to be 
imposed in order to maintain the success of the standard Big-Bang Nucleosynthensis (BBN). Given 
that it does not modify the results in any essential way, other than by moderate factors, we opt to 
keep the analysis as generic as possible, and therefore focus on the thermal production of CDM 
candidates. 

A particle X consists a viable CDM candidate, provided its relic density Clxh 2 can be confined 
in the region [1] 

(a) 0.097 < D. x h 2 < 0.12 for (b) 10 keV <m x < m NLS p (1.2) 

where wjnlsp is the NLSP mass and the lower bound on mx comes from the fact that lower nix's 
cannot explain [19] the observed early reionization [1]. Clearly, the Q-xh 2 calculation crucially 
depends on the assumption on the dominant component of the universe during the decoupling of 
WIMPs or the reproduction of e-WIMPs. In the standard cosmological scenario (SC) we assume 
that primordial inflation is followed by the radiation dominated (RD) epoch. However, our current 
knowledge of the history of the universe before BBN, also allows for other possibilities (see, e.g., 
[20, 21, 22]). Among them, an interesting alternative is provided by the presence of a kination 
dominated (KD) [23] post-inflationary era, which enhances the thermal abundance of WIMPs [20, 
24, 25, 26] and significantly reduces thermal abundance of e-WIMPs [32] with respect to (w.r.t) 
their values in the SC. 

The existence of a KD era is an open possibility in the framework of quintessential scenaria 
(QS). Quintessence [27] (for reviews, see Ref. [28]) is a scalar field, slowing evolving today, which 
can provide the required amount of the present vacuum energy and therefore explain (at least at the 
classical level) the other major component of the universe, the DE. Kination is also an indispensable 
ingredient of the quintessential inflationary scenaria [29, 30, 31]. In recent papers [26, 32] we 
considered the creation of a KD era in the context of the exponential quintessential model [33, 34], 
taking into account a number of relevant phenomenological requirements. Although this model can 
reproduce a viable present-day cosmology in conjunction with the domination of an early KD era, 
for a reasonable region of initial conditions [35], it does not possess a tracker-type solution [24, 36], 
where quintessence is able to reach the required value today starting from a very wide set of initial 
conditions in the remote past. This attractive behavior occurs in models with inverse power-law 
potentials, which are naturally expected in high energy particle physics models [38]. In this way 
the so-called "coincidence" or "why now" problem, related to the fact that the quintessential energy 
density is such that it is dominating the cosmic expansion right now, is addressed. However, within 
the minimal realizations of these models, these positive features do not coexist with the presence 
of an early KD era [36, 39]. 
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In the current work, we reconsider the generation of an early KD era in the context of tracking 
quintessence, switching on a Hubble-induced time dependent mass term for the quintessence field, 
along the lines of Ref. [39] (see also Ref. [40, 41, 42]). The paper is structured as follows: In 
Sec. 2, we show that a KD era within this framework, is characterized by an oscillatory evolution 
of the quintessence field and the barotropic index. Observationally acceptable values for the latter 
at present times can be obtained for relatively law values of the exponent in the potential [43]. 
Other restrictions arising from BBN, the inflationary scale and the DE density parameter can also 
be met in a wide range of the parametric space, which turns out to have a band structure. In Sec. 3, 
we investigate the impact of this KD era on the thermal abundance of WIMPs and e-WIMPs, 
solving both numerically and semi-analytically the relevant Boltzmann equations. We find that, if 
X is a WIMP, Q.xh 2 depends crucially on the hierarchy between the freeze-out temperature and the 
temperature at which the extrema in the evolution of the quintessence field are encountered. On 
the other hand, if X is an e-WIMP, Q.xh 2 is determined mainly at the temperature where the first 
extremum (after the onset of KD era) in the evolution of the quintessence field occurs. In Sec. 4 and 
5 respectively, we study the dependence of Q-xh 2 on the free parameters of the theory, and identify 
the allowed parameter space for WIMPs and e-WIMPs. Our conclusions are summarized in Sec. 6. 
For completeness, we also discuss the cosmology of unstable G within our QS in Appendix A. 

Throughout the text, brackets are used by applying disjunctive correspondence. Natural units 
are assumed for the Planck's and Boltzmann's constants and for the velocity of light (h = c = &b = 
1). The subscript or superscript "0" refers to present-day values (except in the coefficient Vo) and 
log [In] stands for logarithm with basis 10 [e]. Finally, we assume that the domain wall number 
[10] is equal to 1. 

2. Tracking Quintessence 

In this section we outline the several aspects of tracking quintessence (Sec. 2.1) and the various 
observational restrictions that have to be imposed (Sec. 2.2). We then highlight the scalar field 
dynamics in Sec. 2.3 and describe the allowed parameter space of our QS in Sec. 2.4. 

2.1 The Quintessential Set-up 

We present below the equations which govern the evolution of the quintessence field (Sec. 2.1.1) 
and the method we use in order to solve them numerically (Sec. 2.1.2). 

2.1.1 Relevant Equations. We assume the existence of a spatially homogeneous scalar field q (not 
to be confused with the deceleration parameter [3]) which obeys the equation: 

M 4+a b 

q + 3Hq + V q = 0, where V = V a + V h with V a = — — and V b = -H 2 q 2 , (2.1) 

q 2 

is the adopted potential for the field q with M a mass scale and , q [dot] stands for the derivative 
w.r.t q [the cosmic time, t\. The main feature of V a is the existence of tracker solutions, which 
are attractors [36, 37] in field space, while Vb (with b of order unity) expresses a quite generic 
interaction term which arises e.g. due to non-canonical terms of the Kahler potential of q [40, 41]. 
Similar interactions [42] arise due to the thermal effects as well. As shown in Ref. [39] and verified 
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in Sec. 2.4, a mild tuning of the coefficient b enlarges the configuration space that leads to the 
desirable insensitivity to the initial conditions, without modifying the behavior of the field today. 
The Hubble expansion parameter H in Eq. (2.1) is given by 

H = yjPq+PR+pMl V3m P with p q = -q 2 + V, (2.2) 

the energy density of q. The energy density of radiation, pr, can be evaluated as a function of 
temperature, T, whilst the energy density of matter, pm, with reference to its present-day value: 
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PR = -^g P * T 4 and p M R 3 = p M o#o (2.3) 

with R the scale factor of the universe. Assuming no entropy production due to the domination of 
q or another field, the entropy density, s, satisfies the equation 

2tt 2 

sR 3 = s R 3 where s = — g„ T 3 , (2.4) 

where g pi (T) [g s *(T)] is the energy [entropy] effective number of degrees of freedom at temperature 
T. Their precise numerical values are evaluated by using the tables included in public packages 
[46], assuming the particle spectrum of the Minimal SUSY Standard Model. 

2.1 .2 Numerical Integration. The numerical integration of Eq. (2.1) is facilitated by converting 
the time derivatives to derivatives w.r.t the logarithmic time [35] which is defined as a function of 
the redshift z: 

x = \n(R/R ) = -ln(l + z). (2.5) 
Changing the differentiation and introducing the following dimensionless quantities: 

Pr = Pr/PcO, P m = Pm/PcO, V a = VJp c0 , H = H/H Q , and q = q/ V3m P , (2.6) 
Eq. (2.1) turns out to be equivalent to the system of two first-order equations: 

Q = Hq' and HQ' + 3HQ + bH 2 q+ ^H^q 2 + V ' a<g = 0, (2.7) 

where H 2 = ^ _ 1 ^ iy-Q 2 +p R + V a + p M J . (2.8) 

Here, "prime" denotes derivative w.r.t. x and M can be expressed in terms of the dimensionless 
quantities a follows 

M = (( V3m P ) fl t/ Pco) 1/(4+a) with V a = V /q a . (2.9) 
In our numerical calculation, we use the following values: 

p c0 « 8.1 • 10 47 /i 2 GeV 4 , with h = 0.72, p M0 = 0.26 and T = 2.35 • 10" 13 GeV. (2.10) 

We have also H = 2.13 • 10~ 42 /i GeV and from Eq. (2.3) we get p R() = 8.04 • 10~ 5 . 

Eq. (2.8) can be solved numerically if two initial conditions are specified at an initial loga- 
rithmic time ti, which corresponds to a temperature T\ defined as the maximal T after the end of 
primordial inflation, assuming instantaneous reheating. We take g(ti) = 10" 2 throughout our inves- 
tigation, without any loss of generality (see below) and let as a free parameter H\ (which practically 
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coincides with Q(%i)/ V2 since we require a complete domination of kination at early times as we 
describe below). To test our model against observations we extract the density parameters of the 
g-field, radiation and matter, 

= PiliPq +Pr +Pm), where i - q, R and M, (2.11) 
respectively, and the equation-of-state parameter of the g-field, w q , 

w q =y where P q = l -Q 2 -V a - b -H 2 q 2 and p q = ^Q 2 + V a + b -H 2 q 2 , (2.12) 

with P q the pressure of q and P q - P q /p c o- 
2.2 Imposed Requirements 

We impose on our quintessential model the following requirements: 

2.2.1 Constraint of Initial Domination of Kination. As stressed in the introduction, we focus on the 
range of parameters that ensures the initial domination of the ^-kinetic energy. This requirement 
can be quantified as follows: 

Q\ = Q ? (Ti) = 1. (2.13) 

2.2.2 Nucleosynthesis Constraint. The presence of p q has to respect the successful predictions 
of BBN, which commences at about tbbn = -22.5 corresponding to Tbbn = 1 MeV [44]. Taking 
into account the most up-to-date analysis of Ref. [44], we adopt a rather conservative upper bound 
on Q 9 (tbbn)» l ess restrictive than that of Ref. [45]. Namely, we require: 

Qf N = Q 9 (tbbn) < 0.21 (95% c.l.) (2.14) 

where 0.21 corresponds to additional effective neutrinos species 6N V < 1.6 [44]. We do not consider 
extra contribution in the left hand side of eq. (2. 14), due to the energy density of the gravitational 
waves [48] generated during a possible former transition from inflation to KD epoch [29]. The rea- 
son is that inflation could be driven by another field different to q and so, any additional constraint 
arisen from that period would be highly model dependent. Nevertheless, inflation can provide a 
useful constraint for the parameters of our model as we discuss below. 

2.2.3 Inflationary Constraint. Resent data [1] strongly favors that the universe underwent an early 
inflationary phase. Assuming that this phase is responsible for the generation of the power spectrum 
of the curvature scalar P s and tensor P t perturbations, an upper bound on the inflationary potential 
V\ and consequently on Hi can be obtained [49]. More specifically, imposing the conservative 
restriction r = P t /P s < 1, and using the observational normalization of P s [1] we get 

Hi < -^-nipPll 2 ^ //i < 2.65 • 10 14 GeV ^ Hj < 1.72 • 10 56 (2.15) 
V2 



where * means that P s * is measured at the pivot scale k t = 0.002/Mpc. 
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2.2.4 DE-Density and Coincidence Constraint. These two requirements can be addressed if (i) the 
present value of p q , p q o, is compatible with the preferred range of Eq. (1.1b) and (ii) p q has already 
reached the tracking behavior. The two conditions can be implemented [39] if 



where we restrict ourselves to the central experimental value of Q. q o, since, this choice does not 
affect crucially our results on the CDM abundance. 

2.2.5 Acceleration Constraint. A successful quintessential scenario has to account for the present- 
day acceleration of the universe, i.e. [1], 



In our case, we end up with eternal acceleration (w q < -1/3 for x > 0). 

Let us finally note that the results obtained on the age of the universe t and the redshift of the 
transition from deceleration to acceleration, z t , are marginally consistent with the experimental 
data, according to which to = (13.69 ± 0.26) Gyr and z t = 0.46 ± 0.26 at 95% c.l. We do not 
impose the experimental data on these quantities as absolute constraints, due to the observational 
uncertainties in their measurement. 

2.3 The Quintessential Dynamics 

The cosmological evolution of the various quantities involved in the model as a function of x is 
illustrated in Fig. 1 for qi = 0.01, a = 0.5, b = 0.2, T Y = 10 9 GeV (n « -51.2) and Hi = 1.7 • 10 52 
(qBBN _ o.oi and Vq = 9 • 10 8 or M = 4.8 eV). For comparison we also depict by dotted lines the 
evolution of certain quantities as a function of x for b = 0, Vo = 8-8 ■ 10 9 (or M = 8 eV) and the 
same residual parameters. In particular: 

• In Fig. l-(a), we present logp,- versus x for b = 0.2 (solid lines) and b = (dotted lines), 
and for i = q (bold black lines), i = R + M (light gray lines) and i = A (thin black lines). 
Here, p q is computed by inserting in the last equation of Eq. (2.12) the numerical solution of 
Eq. (2.8). The quantity p R+M = p R + p M is given by Eq. (2.3), and p A is the dimensionless 
energy density of the attractor solution (see Sec. 2.3.3 for details). 

• In Fig. l-(b) [Fig. l-(c)], we display q [q'] versus x for b - 0.2 (solid line) and b - (dotted 
line). We observe that for b = 0, q grows to a value greater than mp before it slows down. 
As a consequence, p q overshoots the tracker solution as shown in Fig. l-(a). On the contrary, 
for b = 0.2 this increase of q can be avoided and the tracker solution is reached before the 
present epoch. 

• In Fig. l-(d), we plot - with i - q (black line), R (light gray line) and M (gray line) - 
and w q (dark gray line) versus x. We compute f2, [w q ] applying Eq. (2.11) [Eq. (2.12)]. We 
observe that, in the presence of q, the universe undergoes successively a modified KD era, 
the RD era and then the matter-dominated era until the re-appearance of DE. During this KD, 
w q takes oscillatory values between 1 and - 1 in sharp contrast to the case of a pure KD era 
where w q - 1 - see Ref. [26]. 



(a) Q^o = P q0 - 0.74 and (b) d 2 V{x = 0)/dq 2 - flg, 



(2.16) 



-1.12 < «; 9 (0) < -0.86 (95% c.l.). 



(2.17) 
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x = ln(R/R ) 

(a) 



x = ln(R/R ) 

(b) 
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Figure 1: The cosmological evolution as a function of x = ln(R/R ) for qi = 0.01, a = 0.5, b = 0.2, T\ = 10 9 GeV 
(ti -51.2; and Hi = 1.7 • 10 52 (tf m = 0.01 and V = 9 ■ 10 s ) of the quantities (a) logp, with i = q (solid black line), 
R+M (light gray line) and A (thin line) (b) q (c) q'(= dq/ch) and (d) w q (dark gray line) and Q,/ with i = q (black line), R 
(light gray line) and M (gray line). For comparison we also depict by dotted lines the evolution as a function of x of the 
quantities (a) logp (b) q and (c) q' for a = 0.5, b = 0, Ti = 10 9 GeV (xi -51.2), Hi = 1.7 • 10 52 and V = 8.8 • 10 9 . 



To further facilitate the understanding of the quintessential dynamics we present below a qual- 
itative approach applying the arguments of Ref. [30, 38]. In particular, q undergoes the following 
four phases during its evolution: 

2.3.1 Kination Dominated Phase. During this phase, the evolution of both the universe and q is 
dominated by the kinetic-energy density of q. Consequently, Eq. (2.7) reads: 



Q' +3Q + bHq a and Q = Hq' with H - ^P q ^QI ^2 - bq 2 . (2.18) 
The former equations can be integrated trivially to give: 

q ~ ^sin| V^(x - x T ) + arcsin yj^qi and Q - Q l cos V^(t - x{)e~ %x ~ x ^. (2.19) 

Obviously, for b > 0, q and Q are set in harmonic oscillations during the KD era. Note that for 
b — > 0, we recover the well-known results of a pure KD phase [26] - depicted by dotted lines in 
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Fig. l-(c) and (d). In other words [39], V with b > acquires a (time-dependent) minimum and q 
is prevented from increasing sharply as in the case with b = 0. In particular, q develops extrema at 



(2k + X) U - 



arcsm 



2 qi 



+ ti with k = 0, 1,2,... 



(2.20) 



On the other hand, q' and Q almost vanish for t = x ext (given that q\ ~ 0). Therefore, at x - T ext our 

-2 

approximation in Eq. (2.18) fails instantaneously and pr dominates over Q /2. As a consequence 
the q and Q oscillations become anharmonic and the results of Eq. (2.19) deviate from their numer- 
ical values. These deviations are enhanced as |x - xi| increases and/or H\ decreases. The oscillatory 
behavior of q is the origin of the peaks shown on the curve of p q in Fig. l-(a) and the oscillatory 
form of w q in Fig. l-(b). Nevertheless, the simple formula - see Ref. [26] - which estimates the 
point xkr where the totally KD phase is terminated, is still valid with rather good accuracy, i.e., 



tkr - Tj + 



In yJpql/PRh 



(2.21) 



For b = 0.2 - and the inputs of Fig. 1 - we get numerically [analytically] xkr - -27.5 [xkr - -28] 
which corresponds to 7kr = 0.13 GeV [7kr = 0.2 GeV], whereas for b - the numerical findings 
coincide with the analytic ones (xkr = -28). 

2.3.2 Frozen-Field Dominated Phase. For x > xkr, the universe becomes RD - and so H - 
p R /(l - bq 2 /2) - while the evolution of q continues to be dominated by its kinetic energy density. 
As a consequence, Eq. (2.1) is simplified as follows: 



q" + q +bq+ -qq - 0. 



(2.22) 



Due to the complexity of this equation, it is hard to to obtain a reliable approximate solution. What 
we can say, however, is that, during this period, both q and Q cease to oscillate and freeze at an 
almost constant value. For b - 0.2, logp^ decreases less steeply than for b - and thus, logp^ 
may join the tracker solution in time - see Fig. l-(a). 



2.3.3 Attractor Dominated Phase. In this regime, omitting q and corrections of order b 2 , Eq. (2.1) 
can be simplified as follows: 



q' + b -q- - 7 ^2_ e 3(l + » B )x a where - b = p Boe -3(l + »B)x 



3" 3^ a+1 >p 



(2.23) 



BO 



is the dominant background energy density of the universe with wb - 1 /3 [0] for the RD [matter- 
dominated] era. The solution of (2.23) can be written as 



^ a(a + 2)V \ l/(2+a) 3(Ulll .„ 0; 
/A 1 (9 + (2 + b)a)p M0 



(2.24) 



As a consequence [38], the system in Eq. (2.1) admits a tracking solution since the energy density 
of the attractor 

(2.25) 



p A « p Afe -3(i + «4 p x-M) with fp = a ^zl 



a + 2 
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tracks the dominant p B until t = x A f where the tracking regime finishes and p B - p A . Indeed, p A 
decreases less rapidly w.r.t p B for a > 0, since 

Pa/PB oc e 6(l + » B )x/(« + 2) (2 26) 

As a result, p 9 eventually dominates and the expressions leading to the scaling solution of Eq. (2.25) 
can be neglected. Solving pM = Pa w.r.t. t (since in our cases p B = p M ) we obtain the following 
expression for XAf : 

1 a(2 + a) 1 pmo 

tAf - T<2 m — — - — + -ln^— • (2.27) 

6 (9 + ft (2 + a)) 3 V 

In Fig. l-(a) we depict with a thin solid line the evolution of p A given by Eq. (2.25). For 
the input parameters of Fig. 1 we find that the onset of this phase takes place at x A j - -3.62 and 
terminates at tAf = -0.4 with p Af = 0.88 and ufq -0.81. We check that for x A i < x < XAf, 
p&lp q - 0.96. This fact - in conjunction with the fulfillment of Eq. (2.16b) - ensures a pure 
domination of the attractor solution for a well-defined period, shown in more detail in the subfigure 
of Fig. l-(a). 

2.3.4 Vacuum Dominated Phase. For x > x A f, the universe is dominated by V and so, Eq. (2.1) can 
be written as 

-rf + b -a - ± ~ => a~ £ (l - 0c-™ W+'yP" 9 /PM \ 2fo/9 ) 1/2 28 
9 S 9 3q 9 V^i (9 + b(2 + a)f^\Vo) J " ( ^ 

Using the expression above for q we can estimate w q (0) at present through the formula 

w a (0) « 1 ^— (2.29) 

q 1 + ^'(0) 2 /2 + bq(0) 2 /2 K 

with results ^(0) = 0.42 and u;^(0) ^ -0.88 for the parameters used in Fig. 1. We also obtain 
Zt = 0.76 and to = 13.2 Gyr which are more or less within the experimental limits. 

2.4 The Allowed Parameter Space 

The free parameters of our quintessential model are: 

a, b, M, xi (or T\), qi and H\. 

Agreement with Eq. (2.17) implies < a < 0.6 (compare also with Ref. [43], where a less restric- 
tive upper bound on w q (0) has been imposed). The parameter M can be determined for every a 
through Eq. (2.9), so that Eq. (2.16a) is satisfied. The determination of a and M is independent of 
Xi, qi and H\, provided that the tracking solution is reached in time. This property gives an idea of 
the stability of the tracking solutions. Note that in the case of the exponential potential, studied in 
Ref. [26, 32], any variation on and H\ requires a re-adjustment of Vo so that Eq. (2.16a) is met. 
On the other hand, f2^ BN (which influences the calculation of £l x h 2 ) does depend crucially on Hi 
and Xj (and very weekly on q{). 

In Fig. 2-(a) we illustrate the allowed parameter space of our model in the b - log Hi plane 
for Ti = 10 9 GeV, a = 0.5 and qi - 0.01. In the gray shaded areas, Eqs. (2.13) - (2.17) are 
fulfilled. The upper boundary curves of the allowed bands come from Eq. (2.14). In the overall 
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Figure 2: Allowed (gray shaded) region by Eqs. (2.13)-(2.17) in the (a) b - logHj plane for Ti = 10 9 GeV and (b) 
Ti - log Hi plane for b = 0.2. In both cases we take a = 0.5 and q Y = 0.01 . 

allowed region, we obtain 10~ 6 < H BBN < 0.21. Note, however, that saturation of Eq. (2.14) is 
not possible for 0.08 < b < 0.16. Clearly, the model possesses an allowed parameter space with 
a band structure. If (b, log H\) belongs in a white [gray] band the resulting q after the oscillatory 
phase turns out to be negative [positive] and thus, it cannot [can] serve as quintessence. Let us once 
more fix b - 0.2. For 51.7 < \ogH\ < 53.3, q develops five extrema during its evolution - which 
is of the type shown in the fubfigure of Fig. l-(b) - resulting to q > 0. As log Hi decreases below 
53.3 (where the bound of Eq. (2.14) is saturated), the amplitude of the fifth peak, which appears 
in the ^-evolution (at about t - -24.5) eventually decreases and finally this peak disappears at 
log/?i - 51.7 where the first allowed band terminates. For 48.7 < log/?i < 51.7, q develops four 
extrema during its evolution, resulting to q < 0. As log Hi decreases below 51.7 the amplitude 
of the forth peak which appears in the ^-evolution (at about t - -30) decreases and finally this 
peak disappears at log Hi - 48.7 where the second allowed band commences. Note that in the first 
allowed band f2 BBN increases with Hi but this is not a generic rule (as in the case of a pure KD era). 

Fixing b = 0.2 and letting Ti vary in a range of relatively high values (motivated by the models 
of SUSY hybrid inflation [50]), we depict in Fig. 2-(b) the region allowed by Eqs. (2.13) - (2.17) in 
the Ti - log/Zj plane for a = 0.5 and qi = 0.01. The upper boundary of the allowed region comes 
from Eq. (2.15) for 6.7 < Ti/10 10 GeV < 25.4 whereas for 3.8 < Ti/10 9 GeV < 67 it arises from 
the condition <?(0) > 0. The same applies for the left boundary of the allowed region and several 
parts of its right boundary. On some parts of the latter boundary, Eq. (2. 14) is also saturated. In the 
overall allowed region we obtain 10~ 9 < f2 BBN < 0.21. 

3. Thermal Abundance of CDM Candidates 

We turn now to the calculation of the thermal abundance, O^/j 2 of a CDM candidate, X, which can 
be a WIMP or an e-WIMP. If X is a WIMP we assume that it maintains kinetic and chemical equi- 
librium with plasma, is produced through thermal scatterings and decouples (being non-relativistic) 
during the KD epoch. If X is an e-WIMP, we expect that its relic abundance, due to its early de- 
coupling from the thermal bath, is diluted after inflation at a relatively high energy scale and we 
compute its abundance produced through thermal scatterings and decays during the KD era. In 
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Sec. 3.1 we present the Boltzmann equation that governs the evolution of the X number density 
and then describe the procedure we employ to solve this equation numerically (Sec. 3.2) and semi- 
analytically (Sec. 3.3). 

3.1 The Boltzmann Equation 

Since X's are in kinetic equilibrium with the cosmic fluid, their number density, nx, evolves ac- 
cording to the Boltzmann equation: 

J -((tv) (n\ - nf) if X is WIMP, 

nx + 3Hn x - j c ^ eq2 + ^ ^ T Ki(m . /T) r . if x is e . WTM p j 

where H is given by Eq. (2.2). Let us define the residual symbols of Eq. (3.1) separately: 

3.1.1 The Case of WIMPs. In this case, (crv) is the thermal-averaged cross section of WIMPs 
(hereafter denoted as x's) times the velocity and n x is the equilibrium number density of x, which 
obeys the Maxwell-Boltzmann statistics: 

w * (x) = v^ m * x3/2 e ^' xpi (i) ' where x = h x and Pniz) = 1 + ^jt 1 (3 - 2) 

is a function obtained by expanding the modified Bessel function of the second kind of order n for 
x «: 1 (m x is the mass of x)- Assuming that x's are Majorana fermions, we set g = 2 for their 
number of degrees of freedom. Let us clarify that (crv) can be derived from m x and the residual 
(s)-particle spectrum, once a specific theory has been adopted. Following our strategy in Ref. [26], 
we treat m x and (crv) as unrelated input parameters in order to keep our presentation as general 
as possible. Note, also that far enough from s-poles and thresholds (crv) can be expanded [5] 
non-relativistically as (crv) = a + bx, where a and b are constants. 

3.1 .2 The Case of e-WIMPs. In this case, « eq = £(3)T 3 /n 2 is the equilibrium number density of 
the bosonic relativistic species, m; [</,-] is the mass [number of degrees of freedom] of the particle 
i and K n is the modified Bessel function of the 2nd kind of order n. In the relativistic regime 
(T » mi) C x has been calculated using the Hard Thermal Loop Approximation [12, 14], resulting 
to C x = C*F, where 



r>HT _ 



16f(3)m 



2 

P a =l v -""X 



(3.3) 



l^Ji^i) wifli fa = -4- for Z = a. 

CO) \ 93 32^ /a 



Here, ^ a and M a (with a = 1,2, 3) are the gauge coupling constants and gaugino masses respec- 
tively, associated with the gauge groups C/(1) Y , S C/(2) L and 5 C/(3) c , fe) = (1-634, 1.312, 1.271) 
and (c a ) - (33/5, 27, 72). Note that we include the recently corrected [14] nominator (1.211) in the 
logarithm of C? T . Throughout our analysis we impose universal initial conditions for the gaugino 
masses, M a (M G uj) = M\j% and gauge coupling constant unification, i.e., g a (M G uj) = g G vj with 
^gut - 2 • 10 16 GeV. Eq. (3.3) can be applied self-consistently only for g a (T) < 1 or equivalently 
T > Tq = 10 4 GeV. Towards lower T's, non-relativistic (T «; m,) contributions and X production 
from decays start playing an important role. These contributions have been incorporated in our 
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computation for the case of a, following the formalism of Ref. [32]. In particular, for T «: m,, 
Ca - C^ T has been calculated numerically and r,-'s with i - g, q and B are taken into account 
[13, 15, 32] using the following benchmark values of m,'s: 

(m~ q , m~ g , m s ) = (1, 1.5, 0.3) TeV. (3.4) 

On the other hand, in the case of G, we restrict ourselves to the high temperature regime (no 
formalism for the G production at low T is available to date). We do not include G production from 
sparticle decays in the plasma [16] which can change [9] Q.^h 2 by a factor of about two but does 
not alter our conclusions in any essential way. 

3.2 Numerical Solution 

In order to find a precise numerical solution to our problem, we have to solve Eq. (3.1) together with 
Eq. (2.1). To this end we follow the strategy of Sec. 2.1, introducing the dimensionless quantities: 

nx = n x /p%\ 4 q = r%/p%\ « eq = n*/p%\ W) = V3m P p^(crv), (3.5) 
C x = V3m P p^Cx, rhi = m ; /p^ 4 , T = T/p^ 4 and f , = Ti/H . (3.6) 

In terms of these quantities, Eq. (3.1) takes the following master form for numerical manipulations: 

-_, -_ _ f -(or) (n\ - 4 q2 ) if X is WIMP, 

Hn x + 3Hn x - j ^_ eq2 + ^ M_f ifh 2f Kl ( mi /T) if X is e-WLMP, ( } 

where H is given by Eq. (2.8). The integration of Eq. (3.7) is done until tbbn - -22.5 (an in- 
tegration until gives the same result). In the case of WIMPs, we impose the initial condition 
h x {x x ) = ii x q (x x ), where x x corresponds to the beginning (x - 1) of the Boltzmann suppression of 
n x . In the case of e-WIMPs, we set the initial condition n x (xi) - 0. We use Cx = C% T if Ti » Tq 
and 7kr » Tq. On the other hand, if Ti » Tq and Tkr «: Tq we integrate successively Eq. (3.7) 
from ti to tsusy - - 37 - which corresponds to Tsusy - 1 TeV - with Ca = C? T and then from 
tsuSY to tbbn with C~ a = C^ T . Finally, £lxh 2 is evaluated from the well-known formula: 

&x = Pxo/PcO - mxSoYxo/PcO ^xh 2 - 2.748 • 10 8 Y x0 m x /GeV. (3.8) 

where px = mx nx, Yx - nx/s is the X yield and sq h 2 /p c o = 2.748 • 10 8 /GeV. 

3.3 Semi-Analytical Approach 

To facilitate the understanding of our results we adapt the approach carried out in Ref. [26, 32] 
to our set-up. In particular, re-expressing Eq. (3.1) in terms of the variables Yx, Y^ - n q I s and 
F eq = n eq /s (in order to absorb the dilution term) and converting the derivatives w.r.t t to derivatives 
w.r.t t, we obtain: 

Y^yM^^ 1 -^ if * isWIMP ' where (3.9) 
\g q \CxY^ 2 ifZise-WIMP, v ' 

S ., _ S b _2 . fl + 2 2 /2p R forT«CTKR, 

ifa) - ~1= Wlth s = -777. 9b = 1 - -zq and g q \ 

y/pR p 3 J 4 2 ( 1 for t » t K r. 

In extracting Eq. (3.9) we keep only the first two terms in the parenthesis of Eq. (2.8) and the first 
term of the left hand side of Eq. (3.1) for e-WIMPs (see below). For g\, = 1 (or b = 0) [gb - g q - 1], 
we reproduce the well-known results of the pure KD phase [of the SC]. Let us discuss how we can 
solve Eq. (3.9) in the two cases separately: 
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3.3.1 The Case of WIMPs. An accurate approximate solution to Eq. (3.9) can be achieved, intro- 
ducing the notion of freeze-out temperature, T F = 7\xp) = x F m x (see, e.g., Ref. [26] and references 
therein), which allows us to study Eq. (3.9) in the two extreme regimes: 

• For x <s tp, Y x - Y x q . In this case, it is more convenient to rewrite eq. (3.9) in terms of the 
variable A(x) = Y x (x) - Y x q (x) as follows: 

A' = -Y e / - y <^)A(a + 2F; q ) ^J^ q . (3.10) 

The freeze-out logarithmic time xp can be defined by 

A(x F ) = 6 F Y e x \r F ) => A(x f )(A(x f ) + 2F; q (x F )) - 5f(^f + 2) F; q2 (x F ), (3.11) 

where S F is a constant of order one, determined by comparing the exact numerical solution of 
Eq. (3.9) with the approximate one under consideration. Inserting Eq. (3.11) into Eq. (3.10), 
we obtain the following equation, which can be solved w.r.t xf iteratively: 

( In ^ q )'(x F ) = -y F (^)6 F (6 F + 2)^ q (x F ) y¥/ V^f + D with (3.12) 

* = and ( m y-)V)-'(l-l-| + ^), (3.13) 

where the x - x dependence can be derived from Eq. (2.4). 

• For x » xp, Y x » Y x q and therefore, Y^ - Y^ 2 - Y%. Inserting this into Eq. (3.9) and 
integrating the resulting equation from xp down to 0, we arrive at: 

Y x0 = (y~ F + J F )~ l , where J F = f dx l^- y <o^> and Y xF = (6 F + 1) ^ q (x F ). (3.14) 

Although not crucial, a choice S F = 1 .2 + 0.2 assists us to better approach the precise numer- 
ical solution of Eq. (3.9). 

Inserting Eq. (3.14) into Eq. (3.8) we can obtain Q. x h 2 semianalytically. We verify that this 
result matches well the one found through the purely numerical integration of Eq. (3.7). 

3.3.2 The Case of e-WIMPs. In this case, we focus on the most intriguing possibility, in which 
T\ » ?susy but 7kr «: 7susy- Under this assumption, T\ takes sufficient high values, as suggested 
by the majority of the inflationary modes (see, e.g., Ref. [50]) and f2^ BN naturally takes values close 
to its upper bound in Eq. (2. 14). Yxo can be derived by numerically integrating Eq. (3.9) from x = xi 
until xbbn with Cx = C^F. This is, because Yx is stabilized to its final value, Yxo, not at the onset 
of the RD era - as in the case of a pure KD era [32] - but at a high temperature corresponding 
to x ext for k = in Eq. (2.20). There, the first peak (for x > xi) of the q evolution takes place 
and therefore, /?r dominates over p q instantaneously as explained in Sec. 2.3.1. Consequently, Yxo 
decreases in our QS w.r.t its value in the SC but increases w.r.t its value in the case of a pure KD 
phase. To illustrate this key point we display in Fig. 3-(a) [Fig. 3-(b)] the evolution of log Y x for 
X - G[X - a]asa function of x for T Y - 10 1() GeV [T Y = 5 • 10 7 GeV], Hi - 10 55 [Hi = 9.4 • 10 48 ], 
Mi/2 - 1 TeV [f a - 10 11 GeV] and b = (dashed line) or b = 0.2 (solid line). For b = 0.2 and 
m d = 0.44 GeV [m~ a = 9.5 MeV] we get Q x /i 2 = 0.11 (with X = G [X = a]). 
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Figure 3: The evolution of the logarithm of the X yield, Y x with X = G [X = a] as a function of x (a) f(b)J for 
a = 0.5, T, = 10 10 GeV /Tj = 5 • 10 7 GeVJ, /?, = 10 55 /Hi = 9.4 • 10 48 J, M 1/2 = 1 TeV [f a = 10" GeVJ and b = 
(dashed iine) orb = 0.2 (soiid iine). We obtain Q. x h 2 =0.11 for b = 0.2 and m d = 0.44 GeV fm 5 = 9.5 MeV J. 



In Fig. 3-(a) [Fig. 3-(b)] we observe that for b = 0, where a pure KD era occurs, Yx with 
X = G [X = a] takes its final value, more or less, close to the onset of the RD era for tkr - -28.5 
[tRR - -25.4] - corresponding to Tkr - 0.23 [Tkr - 0.017] -, according to the well-known results 
of Ref. [32]. Namely, for X - G the resulting Yxo (for b - 0) is not so precise, since the used C^ T 
from Eq. (3.3) gives reliable results mainly for T » Tq- On the other hand, for X - a (and b = 0), 
Yx is calculated employing for T < T^vsy, Ca = C^ T obtained from the non-relativistic formalism 
of Ref. [32] with the m;'s indicated in Eq. (3.4). Due to the Boltzmann suppression occurred for 
T < Tsusy, Yx takes its present value at about T - Tsusy (tsusy = -37). On the contrary, for 
b = 0.2 and the inputs of Fig. 3-(a) [Fig. 3-(b)], Yx takes its final value at t - -50 [t - -45] where 
q develops its first - for k = in Eq. (2.20) - extremum. Obviously, Yxo f° r b - 0.2 is much larger 
than its value for b - but still lower than its value in the SC for T - T\. Indeed, for the inputs of 
Fig. 3-(a) [Fig. 3-(b)] we obtain log Y &Q = -5.3 [log Y a0 - -3.94] in the SC. 

Despite the fact that C|! T in Eq. (3.9) has a rather simple form given by Eq. (3.3), it is not 
straightforward to find a general analytical result for the integration of Eq. (3.9). This is mainly due 
to the fact that the integrand includes g q and gi, which depend on the q and Q evolution in a rather 
complicate way. Nevertheless, we can write simple empirical relations which reproduce rather 
accurately (within 5%) our numerical results for X - G [X - a], a = 0.5, b = 0.2, T\ - 10 10 GeV 
or Zi = 10 9 GeV [Tj = 10 9 GeV or Ti = 5 ■ 10 7 GeV]. Namely 

~ ,c 17 T i , \, - / Tj , \](M^/m P m x for X = G, _ 

Cl x h 2 = Uxlog—^+fcix logffi + [a2x\og-\-+b2x) 1 , / , v . (3.15) 
V GeV / V GeV / J [ mpm x /fa for X = a, 

where the numerical coefficients a\x, ct2x, b\x, are given by 

_ / (- 1 - 603 ' 95 - 26 > 1-4223, - 845.78) • 10 11 forX = G, 
(aix, a 2 x, b iX , b 2 x) - | ( _, 250 5, 71.064, 9.5244, - 541.6) • 10 5 for X = a. 

We, thus, conclude that in our QS, for a naturally expected hierarchy among T\, Tq and Tkr, 
the calculation of Q-xb 1 depends exclusively on C^ T and not on C^ T , r,'s and m,'s. 
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4. CDM from Thermal Production of WIMPs 

Employing the formalism developed in the previous section, we can analyze the behavior of £l x h 2 
as a function of the free parameters of the QS (Sec. 4.1) and construct regions allowed by the 
various constraints (Sec. 4.2). 

4.1 The Enhancement of Q x h 2 

The presence of g q > 1 in Eq. (3.12) and, mainly, in Eq. (3.14) reduces /p, thereby increasing 
Q, x h 2 w.r.t its value, Q x h 2 \ sc , in the SC - i.e., setting g q - gj, - 1 in Eq. (3.9). The resulting 
enhancement can be estimated, by defining the quantity: 

a x h 2 - a x h 2 l r 

AQ X = — * • (4.1) 

a x h Isc 



The behavior as a function of our free parameters of AQ X within our QS can be inferred from 
Fig. 4 where we display AQ. X versus m x - see Fig. 4-(a) - for Hi - 6.3 • 10 53 or Hi - see Fig. 4-(b) - 
for m x = 0.5 TeV. We isolate two cases commonly encountered in the analysis of several models, 
fixing (cry) = 10~ 8 GeV 2 (light gray curves) and (crv) - 10~ 10 x GeV 2 (gray curves). We take 
b = (solid lines), & = 0.15 (dashed lines) and b - 0.32 (dotted lines). The chosen Hi results to 
qBBN „ o.Ol, 0.068 or 0.19 for b = 0,0.15 or 0.32 correspondingly, whereas in Fig. 4-(b) H™ N 
ranges from 5 • 10~ 6 to 0.079 for b - or from 0.05 to 0.21 [from 0.016 to 0.21] for b - 0.15 
[b = 0.32]. Due to the band structure of the allowed parameter space of the model - see Fig. 2-(a) 
- not all Q™ N 's are achievable for any b. 

Clearly, for b - we get a pure KD era and our results reduce to those presented in Ref. [26], 
i.e., AQ. X increases when m x or Hi (and consequently Q^ BN ) increases or when (crv) decreases - 
see Fig. 4-(a) and (b). On the contrary, for b + 0, Ml x depends crucially on the hierarchy between 
Tp and text- Given that /p takes its main contribution from g q for x ~ Tp, /p is enhanced - see 
Eq. (3.14) - if tp is lower than T ext and close to it, since g q is suppressed (g q - 1) for x - T ext . 
As a consequence - see Eqs. (3.8) and (3.14) - A£l x diminishes. This argument is highlighted 
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Table 1: The minima of AC1 X , ACl^^, occurring at the freeze-out logarithmic time tp' n or the mass of x m™ m , for 
a = 0.5, T\ = 10 9 GeV, H\ = 6.3 • 10 53 and several (o~v)'s andb's employed in Fig. 4-(a). We also show the range of Tp 
for 0. 1 < m x /TeV < 1 .5 and the logarithmic time T ex t at which the closest to t f peak in the q-evolution takes place. 



by Table 1. There, we list the range of xp for 0.1 < m x /TeV < 1.5 and (crv) = 10 -8 GeV~ 2 or 
(crv) - IO -10 * GeV -2 and the logarithmic time x ext at which the closest to xp's peak in the q- 
evolution takes place for b - 0.15 or b - 0.32 and Hi - 6.3 • 10 53 . Clearly, tp [text] is independent 
of b or H\ [m x or (crv)]. As xp moves closer to x ext , AQ. X decreases with its minimum AQ^| min 
occurring at xp - xj?" 1 or m x = m™ n . Note that xj? 1 " does not coincide with x ext always due to the 
presence of Fp in Eq. (3.14). The appearance of minima can be avoided if xp's happen to remain 
constantly lower than x ext 's. 

Increasing Hi for fixed m x = 0.5 TeV leads to an increase of Q™ N for the parameters of 
Fig. 2-(b). However, the expected increase of AQ X is less effective for (crv) - 10~ 10 x GeV~ 2 than 
for (crv) = 10~ 8 GeV" 2 as shown in Fig. 2-(b). Let us check, e.g., the case for b - 0.32. As Hi 
increases in the first allowed band, shown in Fig. 2-(a), from 52.2 to 53.9, x ext moves form -34 to 
-33.8 and influences AQ. X for (crv) = 10~ 10 x GeV~ 2 more than for (crv) = 10~ 8 GeV -2 . This is, 
because for (crv) - 10" 10 x GeV -2 , we get xp - -33.7 which is closer to x ext 's than xp - -33.1 
which is taken for (crv) - 10~ 8 GeV -2 . Variation of Ti (or equivalently tj) leads to a displacement 
of x ex t's - see Eq. (2.20) - and therefore the minima of AQ. X in Fig. 2-(a) or the limits of the non 
solid lines in Fig. 2-(b) are relocated. However, our results on the behavior of AQ. X remain intact. 

4.2 Restrictions in the m x - (crv) Plane 

Though the post-freeze-out Y x in Eq. (3.14) stays essentially unchanged, residual annihilations 
of ^'s occur up to the present, with several cosmological consequences. Recently, tight upper 
bounds on (crvYs have been reported. These, however, depend on the identity of the products of 
the annihilation of ^'s. To get a feeling of the relevant effects, we adopt the most restrictive bound 
which arises for the annihilation mode of ^'s to e + e~. In particular, the constraints from BBN [51] 
and cosmic microwave background [52] result to the following bounds, respectively: 

TYl 111 

(a) (crv) < 3 • 10" 5 GeV' 2 — and (b) (crv) < 4.4 • 10" 7 GeV -2 — (4.2) 

1 TeV 1 TeV 

Obviously the constraint of Eq. (4.2b) is much more restrictive than this from Eq. (4.2a). To be 
in harmony with the assumptions considered in the derivation of the above bounds, we consider 
hereafter (crvYs independent of T. 

Having fixed the parameters which determine the QS, we can derive restrictions on the param- 
eters involved in the Q. x h 2 calculation. This is done in Fig. 5 where we show the allowed parameter 
space in the m x - (crv) plane for a - 0.5 and Ti = 10 9 GeV. We also take b - and Hi = 6.3 • 10 53 
or 6.2- 10 52 resulting to Q^ BN = 0.01 or 10 -4 respectively in Fig. 5-(a), Hi = 6.3 • 10 53 and b = 0.15 
or Hi = 6.2 • 10 52 and b = 0.32 yielding H™ N = 0.068 or 0.065 respectively, in Fig. 5-(b). 
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by Eq. (1.2) whereas the region above the black solid [dashed] line is ruled out by the ' x = 

upper bound on (crv) from Eq. (4.2a) [Eq. (4.2b)] assuming xx ~ * e+e ■ The conventions 
adopted for the residual lines are also shown. 



In the plots of Fig. 5 the upper bounds from Eq. (4.2a) and Eq. (4.2b) are denoted by black 
solid (left corner of each plot) and dashed lines respectively. On the other hand, the light gray 
shaded regions are confronted with Eq. (1.2). The gray dashed [dotted] lines correspond to the 
upper [lower] bound on Q. x h 2 in Eq. (1.2), whereas the gray solid lines are obtained by fixing Q. x h 2 
to its central value in Eq. (1.2). We observe that Q. x h 2 decreases as (crv) increases. This is due to 
the fact that Q. x h 2 oc 1 /(crv), as can be deduced from Eqs. (3.8) and (3. 14). We also observe that for 
tp far away from T ex t the region allowed by Eq. (1.2) for b ± reaches the one for b - with the 
same H\. However, when tp reaches x ext , Q. x h 2 decreases - as we explain in Sec. 4 - and therefore, 
the required for obtaining Q. x h 2 in the range of Eq. (1.2) (crv) decreases too. As a consequence, 
although the allowed by Eq. (1.2) area in Fig. 5-(b) with b — 0.15 approaches the corresponding 
area in Fig. 5-(a) with the same Hi and violates the bound of Eq. (4.2b) for low rn^'s, it becomes 
compatible with the latter constraint for larger m x s. 

We observe that, due to the dependence of Q. x h 2 on the hierarchy between xp and T ex t for b + 
we can achieve compatibility of the CDM constraint with the bounds of Eq. (4.2b) even for large 
£2 BBN . On the contrary, this can be attained for b = only tuning Q BBN < 10~ 4 . However in 
both cases (b = and b + 0) agreement with the requirement of Eq. (1.2a) implies almost two- 
three orders of magnitude higher (or>)'s than those required in the SC - c.f. Ref. [26]. It is worth 
mentioning that the obtained (crvYs can assist us to interpret [55], through WIMP annihilation in the 
galaxy the reported excess on the positron and/or electron cosmic-ray flux [56], without invoking 
any pole effect [57], ad-hoc boost factor [58] or other astrophysical sources [59]. According to 
preliminary results [60], the best fits to the experimental data can be achieved for the annihilation 
channel^ — > ^lT with m x - (1 - 1.6) TeV and (crv) (1 - 5) ■ 10" 6 GeV" 2 . However, these 
values remain tightly restricted by the CMB constraint [52]. 

5. CDM from Thermal Production of e-WIMPs 

Similarly to the previous section, we analyze the behavior of O-xh 2 as a function of the free parame- 
ters (Sec. 5.1), and we identify the parametric regions allowed by the various constraints (Sec. 5.2). 
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10 10 1 10 10 3 10 10 

m g (GeV) m ff (GeV) 

(a) (b) 

Figure 6: Cl x h 2 as a function of m x (X = G [X = a]) for various Hi 's, indicated in the curves, a = 0.5, b = 0.2 and 
M V1 = 1 TeV [f a = 10 1 1 GeV J (a [b]). We set (a) T Y = 10" GeV (light gray lines) or Ti = 10 10 GeV (gray lines) and (b) 
Ti = 10 9 GeV (light gray lines) or Ti = 5 • 10 7 GeV (gray lines). The CDM bounds ofEq. (1.2) are, also, depicted by 
the two thin lines. 



As emphasized in Sec. 3.3.2, we focus on large TVs which are frequently met in the well-motivated 
models of SUSY inflation (see, e.g., Ref. [50]). 

5.1 £lxh 2 as a Function of the Free Parameters 

By varying the free parameters, useful conclusions can be drawn for the behavior of Q. x h 2 . The 
results are presented in Fig. 6-(a) [Fig. 6-(b)] where we plot £lxh 2 as a function of mx (with X - G 
[X = a]) for a - 0.5, b = 0.2 and M i/2 - 1 TeV [f a = 10 11 GeV], T Y = 10 9 GeV (light gray lines), 
Ti = 10 10 GeV (gray lines) [T Y = 10 9 GeV (light gray lines) and T Y = 5 ■ 10 7 GeV (gray lines)]. 
We use various values for the Hi, indicated in the curves. The CDM bounds of Eq. (1.2) are, also 
depicted by the two thin lines. The results are derived numerically, but they can also be reproduced 
from the empirical expressions in Eq. (3.15). 

For each choice of Ti in Fig. 6, the selected Hi's belong to the allowed region of Fig. 2-(b) and 
yield Q™ N > 10" 5 . In all cases Q.xh 2 is stabilized at T » Tq and therefore, its calculation exclu- 
sively depends on C|! T given by Eq. (3.3). For the TVs and His under consideration, 0™ N and the 
x ext closest to ti, increase with H\. Note that this fact is not explicitly shown in the approximate 
Eq. (2.20). As a consequence, Q.xh 2 (which takes its present value close to this x ex t) decreases as 
Hi increases. This feature is similar to what happens in the pure KD era - c.f. Ref. [32]. Similarly, 
we also observe Q.^b 2 <x 1/m^, whereas f^/i 2 ccm a - see Eq. (3.15). 

5.2 Restrictions in the mx - log H\ Plane 

As already mentioned, X constitutes a good CDM candidate if O-xh 2 satisfies the criterion of 
Eq. (1.2). Enforcing the latter constraint we can derive restrictions in the mx - log//i plane. We 
focus on {T\, HiYs that belong in the first allowed band - depicted in Fig. 2-(b) - of our QS. For the 
sake of comparison, we present results even for b - 0, although the tracking behavior of our QS is 
not attained in this case as explained in Sec. 2.3. 
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5.2.1 Gravitino Cold Dark Matter. We recall that the free parameters in the present case are: 
rag, T\, Mi 12 and H\ with fixed b and a = 0.5. 

In Fig. 7-(a) [Fig. 7-(b)] we plot the allowed (lightly gray shaded) regions in the m^—Hi plane, 
letting Mi/2 vary in the interval (0.5 - 1.5) TeV, for a = 0.5, b = [b = 0.2] and Tj = 10 9 GeV 
or 10 10 GeV. The black solid [dashed] lines correspond to the upper [lower] bound on Cl^h 2 in 
Eq. (1.2a), whereas the gray solid lines have been obtained by fixing D.^h 2 to its central value in 
Eq. (1.2a) for M1/2 = 1 TeV. Having in mind Eq. (3.15), we construct the solid [dashed] line for 
M 1/2 - 0.5 TeV [Mi/2 = 1.5 TeV]. 

The upper [lower] boundary curve (dotted [double dot-dashed] line) of the allowed regions in 
Fig. 7-(a) arises from the saturation of #3 < 1 [Eq. (2.13)]. Recall that #3 < 1 allows employing 
Cjj = C9 T self consistently in our calculation. On the other hand, in Fig. 7-(b), the upper and lower 
boundaries (dashed lines) of the allowed area for T\ - 10 10 GeV arise from the band structure 
of our QS. This is, also, the origin of the lower boundary (dashed line) of the allowed area for 
7i - 10 9 GeV. The upper boundary (thin dotted line) of this area comes from Eq. (2.14). We 
observe that the required m^'s increase with T\ as expected from Eq. (3.15). 

As emphasized in Ref. [32], for b = 0, acceptable Q^h 2, s require a fine tuning of f2 BBN 's to 
very low values. Indeed, in the allowed regions of Fig. 7-(a) for T\ = 10 9 GeV \T\ - 10 10 GeV] we 
have 10" 19 < Q BBN < 10~ 13 [10~ 21 < Q BBN < 10~ 13 ]. Such an unattractive tuning is not needed 
for b - 0.2. In fact, in the allowed regions of Fig. 7-(b) for T\ - 10 9 GeV [T\ - 10 10 GeV] we have 
10~ 4 < Q BBN < 0.21 [10~ 6 < Q BBN < 0.064]. We conclude, therefore, that G is a natural CDM 
candidate within our QS. 
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Figure 8: Allowed (lightly gray shaded) regions in the m- a - log/fj plane for a CDM, 10 10 < / a /GeV < 10 12 , a = 0.5, 
Ti = 10 9 GeV orTi = 5 ■ 10 7 GeV and (a) b = or (b) b = 0.2. For b = we use the m, 's shown in Eq. (3.4). The 
conventions adopted for the various lines are also shown. 

5.2.2 Axino Cold Dark Matter. In considering the candidature of a for the major CDM com- 
ponent of the universe, we have initially to assume that the scalar SUSY partner of a (known as 
saxion) does not decay [13] out-of-equilibrium producing entropy and thereby, diluting Q.ah 2 . We 
then have to ensure the consistency of the hypothesis that a decouples from the thermal bath at a 
temperature Try > T\. To this end, we check that for every T <T\ the following condition is valid: 

H(T) > T- a {T) where T a ~ 6 N F (N 2 - \)q 2 a q\ n e V2. (5. 1) 

Here, is the interaction rate of g's with the thermal bath [54], N F = 12 and N3 = 3 [32]. The 
free parameters in the present case are: m- a , f a , T\ and H\ with fixed b and a = 0.5. 

In Fig. 8-(a) [Fig. 8-(b)] we display the allowed (lightly gray shaded) regions in the m a - 
Hi plane for 10 10 < / a /GeV < 10 12 , a = 0.5, b = [b = 0.2] and T Y = 10 9 GeV or Ti = 
5 • 10 7 GeV. The black solid [dashed] lines correspond to the upper [lower] bound on Og/i 2 in 
Eq. (1.2a), whereas the gray solid lines have been obtained by fixing Q.ah 2 to its central value in 
Eq. (1.2a) for f a - 10 11 GeV. In practice, the solid [dashed] line is constructed for f a - 10 12 GeV 
[fa - 10 10 GeV]. This can be understood taking into account the empirical relations for Qq/z 2 given 
in Ref. [32] [Eq. (3.15)] for b = [b = 0.2]. 

The upper boundary curves (dot-dashed line) of the allowed areas in Fig. 8 come from the 
upper bound on 0™ N in Eq. (2.14). The right boundary (thin line) of the allowed area in Fig. 8-(a) 
arises from the upper bound of Eq. (1.2b) assuming that B is the NLSP, with a mass as in Eq. (3.4). 
Needless to say that this upper bound can be modified if there is another SUSY particle lighter 
than B. The relevant area terminates from below at 7kr - 1 TeV, so that our formulas for C^ T 
in Ref. [32] are fully applicable. The low boundary curves of the allowed areas in Fig. 8-(b) arise 
from the band structure of the parameter space of the QS under study. 
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Contrary to the case of G, we observe that the required m^'s increase when T\ decreases, as 
expected from Eq. (3.15). In the allowed region of Fig. 8-(a) we get 10~ 13 < f2 BBN < 0.21. In 
this case (with b = 0), the £lah 2 calculation is realized employing C^ T corresponding to the m,'s 
indicated in Eq. (3.4). As outlined in Ref. [32] and deduced from Fig. 3-(b), the Qg/j 2 calculation 
in this regime - and therefore, the allowed area of Fig. 8-(a) - is independent of T\, provided T\ > 
Tsusy- On the other hand, in the allowed regions of Fig. 8-(b) for Ti = 10 9 GeV [Ti = 5- 10 7 GeV] 
(and b = 0.2), we have 10~ 4 < Q BBN < 0.21 [0.031 < Q BBN < 0.21] and the Q.- a h 2 calculation 
depends exclusively on C? T , as underlined in Sec. 3.3.2. Larger m a 's are allowed in the case of 
Fig. 8-(a). Obviously, in both cases a turns out to be a natural CDM candidate for a wide range 
of m a 's. However, within our QS (b > 0), this result is insensitive to the low energy s-particle 
spectrum of the theory but depends on T\. 

6. Conclusions 

We studied a quintessential model based on an inverse-power-law potential supplemented with a 
Hubble-induced mass term for the quintessence field, q - see Eq. (2.1). We verified that this term 
ensures the presence of a period dominated by the kinetic energy of q and allows the quintessen- 
tial energy density to develop a tracker behavior sufficiently early, alleviating in this way the co- 
incidence problem. In addition to the numerical treatment of the relevant equations (which is 
mandatory in order to obtain a reliable description of the quintessential dynamics) we presented 
a qualitative but rather comprehensive semi-analytical approach. The parameters of the model 
(a, b, Ti,H{) were confined so as £l q (Ti) = 1 and were constrained by current observational data 
originating from BBN, the present acceleration of the universe, the inflationary scale and the DE 
density parameter. 

We found that < a < 0.6 and that there is a reasonably allowed region in the (b, Hi) plane 
with b mildly tuned to values of order 0. 1 . Extrapolating the results of Ref. [39] to higher tempera- 
tures, we showed that, contrary to the pure KD era, the KD generated in this model is characterized 
by an oscillatory evolution of q and the barotropic index. 

We then examined the impact of this modified KD epoch on the thermal abundance, O-xh 2 , of 
WIMPs and e-WIMPs. Solving the problem numerically and semi-analytically we found that: 

• Q. x h 2 , with^ a WIMP, increases w.r.t its value in the SC. Its increase is not monotonic as 
in the case of a pure KD era, but crucially depends on the hierarchy between the freeze-out 
temperature and the temperature where the evolution of q develops extrema. 

• Clxh 2 with X an e-WIMP (gravitino, G, or axino, a) takes its present value at the closest 
temperature to T\, where q develops its extremum. As a consequence, while Q-xh 2 decreases 
w.r.t its value in SC, it increases w.r.t its value in the pure KD phase, and both G and a arise 
as natural CDM candidates for masses in the range (10~ 4 - 1) GeV. 

We note that, additional BBN bounds might arise due to possible decays of the NLSP Includ- 
ing these effects would have introduced significant model-dependence in our analysis, and would 
be relevant mainly for gravitinos, whose interactions are extremently suppressed. Even in this case 
the changes would be in at the numerical level, while the qualitative features of our discussion 
remain valid. 
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Given that we did not adopt a specific theoretical framework in our approach, we kept for 
simplicity and better definiteness b constant during the various phases of the q evolution - c.f. 
Ref. [39]. However, within supergravity, b may change from phase to phase - see, e.g., Ref. [61] 
-, thus affecting the quintessential dynamics - see Sec. 2.3. We checked that the value of b is to 
remain almost constant (at the level of ±10%) during the RD era which follows the KD era, for 
the tracker solution to be joined in time. On the contrary, if we switch off Vb after the onset of the 
matter domination, our results on Q^o and w q (0) remain more or less intact. Moreover, £lxh 2 , when 
X is a WIMP, increases [decreases] when b decreases [increases] after the KD era - see Eq. (3.14). 
On the other hand, O-xh 2 , if X is an e-WIMP, is not so sensitive to possible alterations of b during 
the post-kination phases, since its magnitude is determined mainly during the KD era. 

Further work [60] is required in order to establish whether the enhancements of Q. x h 2 obtained 
in our QS can explain the reported [56] results on the cosmic-ray fluxes through WIMP annihilation 
in the galaxy. In addition, it would be interesting to check whether the extrema in the evolution of 
q affect the interference between thermal leptogenesis and neutrino masses in conjunction with the 
G constraint - see, e.g., Ref. [62]. 
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Appendix A 

Tracking Quintessence and the G Constraint 

In this appendix, we analyze the implications of our quintessential scheme for the unstable G. In 
this case, G can decay after the onset of BBN, affecting the primordial abundances of the light 
elements in an unacceptable way. In order to avoid spoiling the success of BBN, an upper bound 
on Yjj is to be extracted as a function of and the hadronic branching ratio of G, [63, 64, 65]. 
This is the well-known G constraint. In what follows, we specify some representative values of this 
constraint, taking into account the most up-to-date analysis of Ref. [63]. Note that we here review 
- c.f. Ref. [32] - the m^'s which correspond to different F^TbbnXs, decoding more precisely the 
relevant figures. In particular, if G decays mainly to photon and photino, from Fig. 1 of Ref. [63] 
we can deduce: 



whereas if G decays mainly to gluons and gluinos, from Fig. 2 of Ref. [63] we can deduce: 




and B h = 0.001, 



(A.l) 




(A.2) 
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m G 


min 


Tmax 
KR 




m G 


b = 


0.15 


b = 


0.32 


(TeV) 




(GeV) 




(TeV) 


i?i/10 50 




B h = 0.001 






3.9 • 10 3 


1.8 • 10 4 


1.7- 10 2 


7.9 • 10 3 


0.43 


1.3 • 10 46 


1.1 • 10 5 






%/io~ 17 


0.69 


3.5 • 10 44 


3.9 • 10 6 




0.43 


8.7 


7.6 


22 


15 


10.6 


5 • 10 42 


2.7 • 10 8 




0.69 


4.6 


4 


12 


8.2 






= 1 




10.6 


1.96 


1.7 


5 


3.6 


0.2 


7.9 • 10 46 


2.5 • 10 4 




0.2 


33 


29 


85 


59 


0.67 


1.1 • 10 47 


1.7 • 10 4 




0.67 


4.7 


4 


12 


8.7 


10 


1.5 • 10 44 


8.7 • 10 6 




10 


1.96 


1.7 


5 


3.6 



(a) (b) 

Table 2: The minimum values of Hi, H™'", and the maximum values of 7kr, T^ x , dictated by the G constraint for 
b = and B h = 0.001 or 1 (a) and the obtained Yq(Tbbn) forb = 0. 15 and b = 0.32 and the boundary values of Hi 's of 
the first allowed band depicted in Fig. 2-(a) (b). We use various 's indicated in the tables, a = 0.5,7] = 10 9 GeV and 
Mi/2 = 1 TeV. 



In the SC (where no late-time entropy production is expected) Eqs. (A.l) and (A.2) imply stringent 
upper bounds on T\, for fixed M\j%. For the indicative value M1/2- 1 TeV, we find 

( 10 6 GeV < 0.43 TeV 

Ti < | 2.3 • 10 7 GeV for m d | 0.69 TeV and B h = 0.001, or (A.3) 
1 5.6 -10 8 GeV { 10.6 TeV 

f 1.45 • 10 5 GeV (0.2 TeV 

Ti < J 9 • 10 4 GeV for m d | 0.67 TeV and B h = 1. (A.4) 
1 4.8 • 10 7 GeV I 10 TeV 

Clearly, the upper bound on T\ becomes significantly more restrictive for large B^s and low mg's. 
These restrictions on T\ can be avoided in both the pure kination scenario and our QS. 

Indeed, in the case of a pure KD era the G constraint entails a lower bound on Hi, Hi mm , which 
can be transformed to an upper bound on 7kr, ^™ x , for fixed T\, M\ji and nig - c.f. Ref. [32]. 
Setting Ti - 10 9 GeV and M\ji = 1 TeV, we present in Table 2-(a) the corresponding r?i mm 's 
and T™^ x 's for several m^'s and B^s. Clearly, as the F^'s decrease, the required r?i min 's [r^ x 's] 
increase [decrease]. On the other hand, as we emphasized in Sec. 3.3.2, within our QS is 
stabilized close to the temperature T exp (k = 0) which correspond to x ext (k = 0) - see Eq. (2.20) 
- where the earliest peak of the q evolution occurs. Therefore, we expect that the G constraint 
imposes an upper bound on T exp (k - 0), which is a function of Hi and T\. However, due to the 
band structure of the allowed parameter space of our QS, only certain r?i's are available for each b - 
see Fig. 2-(a). For this reason, we opt to present in Table 2-(b) the F^'s resulting on the boundaries 
of the first allowed band for b - 0.15 and b = 0.32 fixing again Ti = 10 9 GeV and M1/2 = 1 TeV. 
For the selected H{s we obtain T ext {k = 0) a 2.5 • 10 7 GeV [T ext {k - 0) a 8 • 10 7 GeV] for 
b - 0.15 [b = 0.32]. We observe that the obtained F^'s decrease with increasing Hi (as for b - 0) 
and are well below the bounds of Eqs. (A.l) and (A.2) besides for b - 0.32, nig = 0.67 TeV and 
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T, 

Figure 9: Regions in the T\ — \ogHi plane that are simultaneously allowed by the quintessential requirements (gray 
shaded area) - Eqs. (2.13)-(2.17) - and the G constraint for M m = 1 TeV, m d = 0.5 TeV and B h = 0.001 (black lined 
area) or S h = 1 (white lined area). We take a = 0.5, b = 0.2 and q\ = 0.01. 

Bh - 1 where marginally violates the relevant bound. As a consequence, the G constraint can 
be comfortably eluded for the (b, //i)'s used in Fig. 4 and 5 since the employed there H\s belong in 
the ranges of Hi examined in Table 2-(b). Comparing the results of Table 2-(a) and (b) we remark 
that evading the G constraint requires larger His for b ± than for b = - see, e.g., the entries of 
two tables for = 0.67 TeV and Bh - 1 . 

The importance of the KD era in weakening the G constraint within our QS can be also induced 
by Fig. 9, where, in contrast to our previous approach, T\ is variable, whereas is fixed to a 
representative value. Namely, in Fig. 9, we show the regions in the log Ti - log Q™ N plane that 
are allowed by the quintessential requirements - see Fig. 2-(b) -, for = 0.5 TeV, M1/2 = 1 TeV 
and Bh - 0.001 (black lined area) or Bh = 1 (white lined area). We observe that for Bh - 0.001 
the allowed maximal Ti is higher that in the case of Bh = 1. This is because for Bh = 0.001 
we impose F^Tbbn) ^ 1.7 • 10~ 15 , whereas for Bh - 1, we impose Y^(Tbbn) ^ 1.7 • 10~ 16 (in 
accordance with Figs 1 and 2 of Ref. [63]). As a consequence, the maximal allowed T ext (k = 0) 
(4.8- 10 5 -6.8- 10 6 ) GeV for B h = 0.001 is higher than the one (4.6- 10 4 -5.7- 10 5 ) GeV allowed for 
Bh = 1. The same hierarchy holds for 7Ys too. Finally, we remark that the maximal fli's depend 
very weakly on T\. 

Concluding, we can say that although the G constraint is more severe in the present QS than 
in the case of a pure KD era, it remains much weaker than in the case of the SC. As a consequence, 
relatively high values of T\ can be comfortably accommodated in both former cases. 
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